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Abstract 

We propose a non-equilibrium continuum dynamical model for the collective 
motion of large groups of biological organisms (e.g., flocks of birds, slime 
molds, etc.) Our model becomes highly non-trivial, and different from the 
equilibrium model, for d < d c = 4; nonetheless, we are able to determine 
its scaling exponents exactly in d = 2, and show that, unlike equilibrium 
systems, our model exhibits a broken continuous symmetry even in d = 2. 
Our model describes a large universality class of microscopic rules, including 
those recently simulated by Viscek et. al. 
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The dynamics of "flocking" behavior among living things, such as birds, slime molds 
and bacteria has long been a mystery. Recently, a number of simple numerical models 
that exhibit such behavior have been studied [0,|TJ. For example, reference || considers a 
synchronous, discrete time step rule in which an individual "bird" in a group of "birds" 
determines its next direction of motion on each time step by averaging the direction of 
its neighbors in certain area, and then adding some zero mean noise, while keeping the 
magnitude of its velocity constant. Their simulations in two dimensions find a transition 
between an ordered phase in which the mean velocity of the flock < v >^ and a disordered 
phase with < v >= as the strength of the noise is increased. 

The above two dimensional model is very similar to the 2D XY model [Q,fp because the 
velocity of the "bird", like the local spin of the classical XY model, also has fixed length 
and continuous rotational symmetry. Indeed, it is easy to see that, in the limit that the 
magnitude of the velocity goes to zero, on each time step the "birds" are just picking a new 
direction, but never actually move, the model reduces precisely to the Monte Carlo dynamics 
of a two dimensional XY model, with the (small) bird velocity playing the role of the XY 
spin. Since the 2D XY model does not exhibit a long range ordered phase at temperatures 
T > (due to spin wave fluctuations), the long range ordered state observed in reference 
seems very surprising. Indeed, in light of the Mermin- Wagner theorem [f| for equilibrium 
systems, its existence must depend on fundamentally dynamical, non-equilibrium aspects of 
the model. In this paper, we show, using a continuum dynamical equation which describes 
a large universality class of related dynamical models, that this is indeed the case. In 
particular, we explicitly demonstrate: 

1) that our model differs from the equilibrium system for spatial dimensions d < 4, 

2) We can calculate the scaling exponents of this model exactly for d = 2, and 

3) the model does, indeed, have a stable spontaneous symmetry broken state even in two 
dimensions. 

Our starting point is the continuum equations of motion (EOM): 
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d t v+ (v- V)u = av-(3\v\ 2 v- VP + D L V{V -v) + D 1 V 2 v + D 2 (v- Vfv + f (1) 

|^ + V»)=0 (2) 

where f3, Di, D 2 and Dl are all positive, and a < in the disordered phase and a > in 
the ordered state. The left hand side of eq. (1) is just the usual convective derivative of the 
coarse-grained velocity field v. The a and (5 terms simply make the local v have a non-zero 
magnitude (= \Jaj ~/3) in the ordered phase. D L 1)2 are diffusion constants. The Gaussian 
random noise / has correlations: 

< fiftWjF,*) >= A5 l3 5 d (f- 7)5(t - t') 

where A is a constant, and i , j denote Cartesian components. Finally, the pressure 

oo 

P = P(p) = Y / Vn(p-PoT, 
n=l 

where po is the mean of the local number density p(r). The final equation (2) reflects 
conservation of birds. 

The essential difference between our model and the equilibrium XY model is the existence 
of the convective term in our model, which makes the dynamics non-potential and further 
stabilizes the ordered phase. A heuristic argument for the stabilizing effect of the convective 
term can be given if we consider our model in Lagrangian coordinates. In those coordinates, 
the convective term drops out, and the interaction between the velocity field is local at 
each instance. However, at different times, the "neighbors" of one particular "bird" will 
be different depending on the velocity field itself. Therefore, two originally distant "birds" 
can interact with each other at some later time. It is exactly this time dependent variable 
ranged interaction which stabilizes the ordered phase. 

To treat the problem analytically, it is more convenient to use the Eulerian coordinates 
as in eqs. (1,2). In the rest of our paper, we concentrate on studying the symmetry broken 
phase, where a > 0. We can write the velocity field as v = voX\\ + 5v, where voX\\ =< v > 
is the spontaneous average value of v We will ignore fluctuations in the magnitude \v\ from 
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its optimal value of \J~a/P since they decay in a finite time (of order ^). Choosing our units 



of velocity so that ya/(3 = 1, we can now write the velocity as: v = (v±, y 1 — \v±\ 2 ) ~ 
(v±, 1 — ^ | iT_L | 2 ) , provided \v±\ 2 << 1. 

Shifting to a co-moving coordinate frame moving with velocity voX\\, v = (v±, — ||i>j_| 2 ), 
and the convective term becomes: (v± ■ V±)v± — ||#±| 2 djji>j_. We will neglect the second 
term, this will be justified a posteriori. The equation of motion then becomes: 

d t v± + X(v± ■ V±)v± = -V±P + D L V ± (V± ■ v±) + DiVlvx + D {l df { v ± + f ± (3) 

^ + PoV ± • v± + AV • (vdp) = (4) 
at 

with D\\ = Di + D 2 and 5p = p — po- The book-keeping coefficient A = 1 in the physical 
case. 

We first study the linearized EOM. We rescale lengths, time, and the fields v± and Sp 
according to 

x± — > bx±; x\\ — > b^x\\; t — > b z t; v± — > 6 x -u_j_; 5p — > 6 Xp 5p (5) 

we choose the scaling exponents to keep the diffusion constants D\\, D± = D\ + D L 
and the strength A of the noise fixed. The reason for choosing to keep these particular 
parameters fixed rather than, e.g., a±, is that these four parameters completely determine 
the size of the equal time fluctuations in the linearized theory, as can be seen by solving 
that theory in Fourier space: 

< v i {Q,t)Vj {-q,t) >= — [jjr — 2 — j ' 2 . 2 + , 2 - ' 2 ) 
2 {D ± q{ + D\\q^)q{ (D iq 2 ± + D\\qfi)ql 

The exponents for the linear theory can be determined very easily: z = 2, ( = 1, x — l — d/2, 
and Xp — X because the density fluctuations of 5p are comparable to those of v±. Therefore, 
the linearized theory implies that v± fluctuations grow without bound (like L x ) as L — > oo 
for d < 2, where the above expression for x becomes positive. This implies the loss of long 
range order in d < 2. 



Making the rescalings as described above, the other parameters in the model scale as: 
A ~ 6 7 *A, a n ~ W n a n with 7 A = % + 1 = 2 - d / 2 and 7„ = z - \ + n\ = n + (1 - The 
first of these scaling exponents to become positive with decreasing d are 7a and 72, which 
both do so for d < 4, indicating that the X(v± ■ V)v± and 02V_i_(5p 2 ) non-linearities are 
both relevant perturbations for d < 4. So for (i < 4, the linearized hydrodynamics will break 
down. 

An e = A — d expansion will obviously not be of much use in our problem in d — 2. 
But fortunately, because of the various symmetries in eqs. (3,4), we can obtain the exact 
scaling exponents in d — 2. First of all, the reduced equations of motion eqs. (3,4) have a 
"Galilean Invariance" ||: i.e., if we let: v±(r,t) — > v±(r,t) + v±fi and simultaneously boost 
the coordinate: x± — > x± —Xv±flt, the equations (3,4) remains invariant for arbitrary values 
of v±fi. This implies that there are no "graphical" renormalization of the nonlinear vertex 
(v± ■ V±)v±, it can only renormalize by rescaling. Furthermore, in precisely two dimensions, 
D\\ and A are also only renormalized by rescaling. To see this, note that in two dimensions 
v± has only one component (call it v x ), which can be written as v x = d x h. The equations of 
motion (3,4) can then be rewritten in terms of h: 

d t h + ^\V ± h\ 2 = -V ± P + D ± V 2 ± h + D {l d 2 h + rj (6) 

^ + poVi/i + AV ± ■ (5pV ± h) = (7) 
where the new effective noise rj = has long ranged correlations. In Fourier space: 

< V (q,uj) V {-q", -a/) >= A W % V 

Q± 

This model only corresponds to our original model (3,4), and hence only describes "birds", 
in d = 2. However, we will analyze this model (6,7) in arbitrary spatial dimensions d, with 
the goal of understanding it in the physical case d = 2. 

It is easy to see that Dn and A cannot be renormalized in this model. A cannot be 
renormalized because it is the coefficient of a non-analytic, long-ranged noise-noise correla- 
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tion function. The non-linear couplings A and o~ n in (6,7), being analytic (i.e., local in space 
and time), can therefore not generate such non-analytic correlations. 

The diffusion constant Dm likewise cannot be renormalized, because any such renormal- 
ization must clearly involve at least one A with at least one external h leg, and hence at 
least one power of q±. This can (and does) renormalize D±, but cannot renormalize D». 

Implementing the dynamical renormalization group |)||| we find, quite generally: 

dDU (z-2C)D {l 

(* + l-d-C-2x)A 
( X + z-l)X 
(z-2 + G ± ({g m }))D ± 
(z-l + G p ({g m })T) Po 



dl 
dA 

Hi 

dX 

m 

dD ± 

~dT 

dpo 

dl 
do\ 

~dT 

d(J n , . f ^ 1 . G n ({g m }) 

— = (z + {n- l)x-l + )a n (8) 

dl g n 



l + G 1 {{g m })T)a 1 



with the parameter T = ^-L_ , and the effective nonlinear coupling constants g\ = D 5/t^i/4 

n-l n ' 

and g n > 2 = Z+f^n-r —^r ■ The G± p n 's denote the non-vanishing graphical corrections to 
D±, p , and a n , respectively. Note that they explicitly depend only on the coupling constants 
<7m's, by construction. Note that the absence of graphical corrections to Du, A, and A is 
exact to all orders in perturbation theory, as discussed earlier. Since we seek a fixed point at 
which A, A, and D\\ remain fixed, we get the following three exact constraints on the three 
exponents: 

x + z = l;z = 2C;d + C + 2x = z + l (9) 

whose solution is: 

2(d+l) . d+1 3 -2d 

Z = ^^> C = ^ ! X= — (10) 

If we combine the above RG eqs.(8) , we can obtain the RG flow equations for the 
effective coupling constants g n , and the parameter T: 
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f = -(1 - G ± ({g m }))T - ^(G 1 ({ 5m }) + G p {{g m }))T 2 (11) 

^ = ^(4-d-§G ± ({( 7m }))»i (12) 

= ^(2n + (1 - n)%„ + G„({«U) + - GOr - (13) 

from which we see that g 1 and g 2 become relevant below d — 4, while all g ra>2 are irrelevant 
near d = A. Hence we can neglect all g„>2's, and, therefore, all of the u n >2 , s as well, at least 
near d=4. The parameter T only becomes relevant below d = 1.5 where x > 0, i.e., when the 
ordered phase disappears. We have calculated the graphical corrections G± and G 2 to one 
loop order near d = 4, and obtain: G± = + #2) ^, G2 = 0. We do not know whether 

the vanishing of G2 to this order is the result of some symmetry of the problem that we have 
failed to recognize, in which case G2 would vanish to all orders, or if it is purely coincidental. 
In either case, to one loop order, inserting these results for the graphical corrections into the 
recursion relations (16) yields a fixed line (actually a fixed hyperbola): (y + g2)g\ = 
This summarizes our picture for model (6, 7) in 4 — e dimensions. What happens as we 
move down down to two dimensions, which is the only dimension in which the model (6,7) 
actually describes "birds"? If G2 vanishes to all orders in e, we will still get a fixed line, 
as in 4 — e dimensions, all the way down to d — 2, although its position will be shifted 
(and it might become curved) away from that given by the one loop calculation. If G 2 does 
not remain zero, then the fixed line collapses to a fixed point. In either case, the scaling 
exponents continue to be given by equation (10), since those results depended only on the 
symmetries of the model. 

So in d = 2, the exponents are given by equation (10), i.e., z — |, C — |) an d X — ~ |- 
These exponents can be checked experimentally (or from simulations) by measuring, e.g., 
the density-density correlation function, which is given, in Fourier space, by: 

< \p(q,uj)\ 2 >= (cj2 _ c2g 2 )2 + u2(£>*(£ oj)ql + D^f (M) 



where c = ^0\p is the speed of sound, and D± is the renormalized diffusion constant. As 
a function of u, this correlation function (like all of the correlation and response functions 
for this problem) has two sharp peaks at u — ±cg_i_, of width D^q 2 ^ + -D||<?j]. Thus, in the 
frequency regime containing most of the weight of the correlation function, D^(q±,qu,u) 
can be evaluated at uj = cq±. Using standard renormalization group arguments and the 
recursion relation (8) for Dj_, we find that: 

Df(q ± ,q lb u ^ cq ± ;X, Po ,a n ) = q z ± - 2 f(^) (15) 

Similar RG arguments yield the finite size scaling of the real-space, real-time rms fluctuations 
of v±: 

L 2 L 

< \v±(f,t)\ 2 >= constant — L 2 _^g(—^-) = constant — L ± 5 g( — ^) (16) 

L ± " L\ 

where L± and Ln are the spatial dimensions of the flock perpendicular to and along the 

mean direction of motion, respectively, and in the last equality we have used the value of \ 

in d — 2. Since this goes to a finite constant as L —>■ oo, we see that long ranged order is 

stable in this model in d = 2, as we claimed in the introduction. 

Now that we have obtained all the exponents, we need to return to our original model 

(1,2) and verify a posteriori all of our assumptions. In particular, we must show that it 

was valid to neglect |'5j_| 2 9||'yj_, which, under the rescalings eq.(3) scales like ~ b 2x+z ~ i ' = b s . 

Using the linearized results for x, % and (, we get: 5 = 2\ + z — £ = 3 — d which is clearly less 

than zero near d = 4. Does it remain < down to d = 2? Experience with, e.g., the 4 — e 

expansion for the 4 theory of critical phenomena suggests that it does. In that problem, a 

6 perturbation also has linearized RG eigenvalue 3 — d. This term nonetheless appears to 

remain irrelevant all the way down to d = 2, judging by the success of extrapolations of 

4 — e results for the Ising model down to d — 2. The apparent contradiction between this 

result and the eigenvalue 3 — d, which of course becomes positive for d < 3, is that graphical 

corrections of 0(e) to this result occur, and keep the eigenvalue negative down to d — 2. 

It seems just as safe to assume that this happens here as in 4 theory, and so we strongly 

suspect that it does, and that our results for the exponents do hold exactly in d = 2. 
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Even in the wildly unlikely event that the cubic vertex does become relevant above d = 2, 
however, we can still show that our model has long-ranged order in d = 2. If the cubic vertex 
does become relevant, we can no longer obtain the exact scaling exponents in d = 2, because 
both A and D\\ are now renormalized. 

However, not all the scaling relations are lost. The random force is still unrenormalized 
since even the contribution from the new vertex |£j_| 2 d||?7]_ are proportional to gjj, which 
still vanishes as \q\ — > 0. Furthermore there is a new scaling relation coming from the 
rotational invariance, i. e., since the direction in which we choose to break the symmetry of 
v was arbitrary, we must, even after renormalization, be able to resum the nonlinear terms 
(v± ■ V)v± and | , uj_| 2 9||'Ui vertices into the form (v ■ V)v. This requires that the graphical 
corrections to (v± ■ V)v± and |#i_| 2 9i|?;j_ be the same. To find a fixed point, therefore, we 
must have their rescalings to be the same as well. This leads to a new scaling relation: 
2x — 1 = 3% — C, or x — C — 1- Taken together with the last exponents relation in (9), this 
leads to the following scaling relation between z and X'- X = Now we expect on physical 
grounds that z < 2 for all dimensions d < 4, since physically, the motion of the "birds" 
enhances the mixing, and we know the corrections due to this effect diverge below d — 4. 
Hence we expect hyperdiffusive behavior, which implies z < 2. Then the scaling relation 
X = implies x < 0, i.e., true long range order in d = 2. 

Numerical simulations [[| indeed find a long range ordered state in the low "temperature" 
regime, in agreement with our predictions above. Detailed study of the correlation functions 
to test our predictions for the scaling exponents (e.g., measurements of the p — p correlation 
function in eqs.(14)) would clearly be of great interest. 

Considerable work remains to be done on this model. The properties of the low tem- 
perature phase of our original model (1,2) in d > 2 remain to be determined. Since the 
symmetries which prevent the renormalization of the noise strength A and the diffusion 
constant D\\ are lost in d > 2, it is no longer possible to obtain exact exponents. However, 
an e expansion on the full model (1,2) should give quite accurate exponents in d = 3. 

There is also the question of the transition from the ordered to the disordered state. 
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Without the convective vertex, our model (1,2) is just model A dynamics for a 4 theory, as 
studied by Halperin et al. However, we can show that, as in the low temperature phase, the 
convective vertex becomes relevant at the transition in d = 4 as well. An e expansion study 
of this problem is also currently underway We will include these subjects and the detailed 
account of this letter in future publication || . 

We are grateful to T. Viscek for introducing us to this problem, and providing us with 
an early draft of reference 2. We also thank P. Weichman, G. Grinstein, and D. Rokhsar for 
many valuable discussions. 
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